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Abstract. - We explore an efficient way to numerically evaluate the response of the surface 
stress of a metal to changes in its superficial charge density by analysis of the strain-dependence 
of the work function of the uncharged surface. As an application we consider Au(lll), (110) and 
(100) surfaces, employing density functional calculations. The sign of the calculated response 
parameter can be rationalized with the dependence of the surface dipole and the Fermi energy on 
strain. The numerical value falls within the range indicated by experiment. The magnitude can 
explain the experimentally observed volume changes of nanoporous materials upon charging. 



Recent experiments reveal a macroscopic expansion or contraction when high surface area 
metals with nanometer-sized porosity are electrically charged, the signature of changes in 
the surface bond forces when a space-charge layer is generated [1,2]. The parameter which 
quantifies these forces in the continuum description of solid surfaces, the surface stress, /, 
is a topic of current interest in surface science, since it is intricately related to the surface 
electronic structure and bonding, and since it is relevant for reconstruction as well as for the 
stress in thin film devices [3,4]. Furthermore, the impact of the energy-conjugate quantity to 
/, the tangential strain e, on the surface electronic structure is being recognized as central for 
the catalytic activity of metal surfaces [5,6]. The response of / to changes in the superficial 
excess charge density (per area), q, relates to fundamental issues in electrochemistry, such as 
electron transfer in surface-adsorbate bonds or microscopic processes in the electrochemical 
double layer [3,7-9]. Experiments finding a stronger response when there is less adsorption 
point towards the effect of q on the bonding in the metal surface as a decisive factor [10, 
11]. Yet, the microscopic processes linking q to forces and relaxation at metal surfaces are 
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poorly understood. A recent density functional theory (DFT) simulation of Au cluster ions 
highlights the role of stretch (relaxation of the surface in normal direction) in response to 
charge [12]. Qualitatively, the relaxation might be rationalized as driven by out-of-plane 
Hellman-Feynman forces between the excess charge and the surface layer of atoms. It was 
argued that the outward relaxation feeds back into the (in-plane) surface stress owing to the 
transverse contraction tendency of solids [12]. This hypothesis, however, awaits verification 
by experiments or simulation using single crystal surfaces of different orientation. 

Here, we go a first step into this direction and explore the response of / to g for surfaces 
of gold. We begin by outlining an ab initio approach, which allows to investigate the response 
with good accuracy. In close analogy to a procedure which has recently been proposed for 
experimental work [13], we exploit a Maxwell relation that permits studying the response 
without explicit introduction of an excess charge into the calculation. For a general discussion 
of Maxwell relations with an eye on applications to interfaces, we refer the reader to the text- 
book by Ibach, Ref. [14]. Here, we apply the concept to gold surfaces using DFT calculations, 
and we give a brief discussion of the relevant microscopic processes. 

As a model situation, we consider a metal plate extended laterally and of finite thickness. 
The energetics of the extended bulk metal and the excess at the surface are represented by 
energy densities, 5" and ip, respectively. In experimental studies and in continuum theory 
the densities are commonly measured per volume or per area, whereas ab initio studies often 
specify energies per atom or per crystallographic unit cell. Both conventions may be treated 
by a common set of equations, provided that the surface area. A, and plate thickness, d, 
of the continuum picture are measured in coordinates of the undeformed lattice (Lagrange 
coordinates; see Ref. [15] and references therein), a convention which we adopt here. More 
specifically, we take ^ and -0 to obey 

Edie,q)^Ad^{e) + AMe,q) (1) 

where E is the plate's total energy and the subscript - which we shall drop henceforth for 
brevity - accounts for a possible dependency of effective surface properties on the number of 
layers in thin plates. We allow for the plate to be elastically strained in the tangent plane, 
while it is free to relax along its normal, with no external forces acting on its surface. The 
tangential strain, e, represents the relative change in surface area - as measured in laboratory 
coordinates - owing to the in-plane elastic deformation. Our particular interest is in the 
dependence of V' on g and on e. A scalar surface stress is defined via 

f{e,q) = di;/de\g. (2) 

Generally, the forces by which the surface interacts with the underlying bulk crystal are 
anisotropic. For planar surfaces they can be described by a 2 x 2 tensor, of which / is half 
the trace. For reference below we note that, by symmetry, the surface stress i s iso tropic on 
the (111) and (100) surface of an fee lattice, but it may be anisotropic on (110. If^)! 
The response of / to charging defines a scalar surface-stress-charge coefficient [15], 
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(^) To include the anisotropy, let ^ and i/i depend on the tangential strain tensor, e, instead of e [15]. The 
surface stress is here obtained as a second rank tensor in the plane, s(e,g) = dip/de\q. The stress-charge 
coefficient of Eq. JSjl is then also a second rank tensor, <; = ds/dq\^ . The appropriate form of the Maxwell 
relation, Eq. JSJ, is ds/dq\e = —qo~^dfi/de\q. 
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Table I - Properties of unreconstructed Au-surfaces at e — as obtained in this work: free energy, 
surface stress, f, work function, W, and surface stress response, and f have been determined 
via Eqs. flp and fH), extrapolating 1/d — > 0. Published values obtained using LDA are shown for 
comparison, and are marked by reference. 





[J/m"] ([eV/cell]) 


/ [J/m^] 


W [eV] 


<>-[v] 


(Ill) 


1.114 (0.502) 


3.317 


5.65 


-1.86 




1.25 [24] 1.39 [25-26] 


2.77 [24] 


5.63 [29] 






1.68 [27] 


2.56 [27] 






(100) 


1.343 (0.700) 


2.723 


5.56 


-0.90 




1.62 [25-26] 1,44 [28] 


3.11 [28] 


5.61 [28], 5.53 [29] 




(110) 


1.372 (1.010) 


2.020 


5.42 


« 




1.75 [25-26] 




5.41 [29] 





Let us introduce the electron chemical potential, /i = dE/dn\e with n = —Aq/qo the number 
of electrons in the system, where go denotes the elementary charge, (7o=l-6022 x 10^^^ C. By 
inspection of Eq. ^ it is seen that 

/i = -qodip/dqle- (4) 

The Maxwell relation for ip{e,q), 

df/dq\e = -qa^'^^^i/^e\q, (5) 

relates <j to a change in /i when the surface is strained tangentially at constant q. 
The work function, W , is related to /i by 

A* = -W = - £f, (6) 

where 4>oo is the asymptotic vacuum value of the electrostatic potential and is the energy of 
the highest occupied single particle state (Fermi energy) of the metal plate [16]. Therefore, as 
an immediate consequence of Eq. ([5]), the linear response quantity d f / dq\q=a can be obtained 
by calculating the work function change when the uncharged slab is strained uniformly in the 
plane, and using Eqs. (O and ([6]). This is a routine calculation with standard codes. The 
advantages are twofold: a) Extending codes to include excess surface charges is not required 
at this point, b) Directly determining / using the definition Eq. ^ comes with technical 
difficulties: in order to obtain ip from total energy computations, total energies from slab- 
calculations with different thicknesses are to be compared. Taking finite energy differences 
between cells with different sizes keeping very high numerical accuracy (needed for derivatives) 
can be computationally expensive, because the number of basis functions varies with the cell 
size. When using Eq. ((5]), explicit knowledge of the free energy is not required. 

Next, we apply our method and determine <, for Au surfaces (111), (100) and (110) from 
Eq. ([5]). To this end we performed DFT calculations within LDA [17] using a mixed-basis 
pseudopotential (MBPP) approach [18-21]. The core- valence interactions in Au were de- 
scribed by norm-conserving ionic pseudopotentials, and the valence states were represented 
by a mixed basis of plane waves {Epyj — 20 Rydberg) and five local d-orbitals per Au atom (r;o 
= 2.5 Bohr). Monkhorst-Pack k-point meshes were used for Brillouin-zone integrations [22]. 
Our present computational LDA-MBPP approach is the same as for earlier ab-initio studies 
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Fig. 1 - (color online) (a)-(c) Electron chemical potential, /i,of the Au(lll), (100) and (110) slab 
models at various tangential strains, e, as a function of the number of atomic layers A^. (d) as a 
function of e for the thickest samples in each orientation. Lines: fit by linear (111, 100) and quadratic 
(110) functions. 



of gold surfaces by Ref. [23] (and references therein). We consider a slab model with periodic 
boundary conditions in all dimensions, containing two surfaces exposed to vacuum with an 
odd number, TV, of atomic layers (5, 7, 9, ■ • ■) and a vacuum region thicker than 13 A, which 
corresponds to at least 6 atomic layer thicknesses. (We have concentrated on the case N odd, 
to avoid even-odd effects, which impair the extrapolation N — > oo.) The thick vacuum region 
is introduced to prevent undesirable interactions of the surfaces via the periodic boundary 
conditions. As a reference configuration for the strain we take the film with lattice parameter, 
a, identical to the computed equilibrium value of bulk Au, uq = 4.085A. The isotropic in-plane 
strain e is applied (i.e., equal strain components of magnitude e/2 along two orthogonal direc- 
tions in the plane, giving a relative area change of magnitude e), and the layers are relaxed, 
keeping their in-plane dimensions fixed while allowing free relaxation out-of-plane, until all 
the relevant atomic forces are below 0.01 eV/A. 

The DFT calculations directly give the value of the highest occupied electronic level. 
The electrostatic potential in the vacuum region saturates rapidly with increasing distance to 
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the metal surface. Thus, the value for the electrostatic potential at the center of the vacuum 
layer gives an excellent estimate for (/loo- Then, W and fi follow from Eq. Table 1 

compares published LDA-data for ^, W, and / for surfaces of different orientation to our 
values. It is seen, that the trends agree well. 

Figure [Ija) shows /i as a function of N for (lll)-oriented slabs for various e. The most 
obvious observation is that varies systematically, and essentially linearly, with e. Fig. [TJa) 
also reveals an interesting detail: the work function (/i — —W) at nonzero strain is very 
sensitive to the layer thickness in thin slabs. In particular when N is doubled from 5 to 11 
at e = 0.04, the resulting shift can be as important, as the change induced by doubling the 
strain. The data converge reasonably well for iV>9 atomic layers (20 A). For the k-point mesh 
density, convergence is achieved with 24x24x2 mesh points. 

The values of /i for the thickest film, A^=13, are plotted over the strain, e, in Fig. [Hd). 
The linear variation is confirmed, and the slope is determined as ^ = — 1.86±0.016V. The sign 
is negative, indicating that the slab tends to expand laterally when positively charged. (Recall 
that a positive / represents a tendency of the surface to laterally compress the underlying 
crystal.) 

Our evaluation of the Au(lOO) surface follows exactly the same lines and /i has been plotted 
in Fig. [IKb). Again, there is a linear behavior, and we obtain ~ —0.90 ± 0.02V, a factor of 
two smaller than for Au(lll). 

In contrast to the previous cases the electron chemical potential of the (110) surface is 
very difficult to obtain. The thick layer limit is not reached, even at A^=17, as shown in Fig. 
[TJc). The strain-dependence of fi for the (110) films differs noticeably from that of the (111) 
and (100) films. Taken at face value, the data (Fig. [Hd)) imply a quadratic, rather than a 
linear variation. 

The most important systematical source of error in our calculation is due to the approxi- 
mation in the exchange correlation functional used for DFT. We have used LDA here, which 
is a controlled approximation for the homogenous electron gas, i.e. in the absence of density 
inhomogeneities. It is well known, that the application of LDA to vacuum surfaces is not 
controlled, in the sense that the expansion parameter (density gradient times Fermi wave- 
length) is no longer small. Therefore it is not clear that gradient terms leave unaffected the 
quantitative results presented in Table [H Further calculations testing functional dependencies 
are necessary and our future work is devoted to this issue. 

Despite of this disclaimer, we believe that, qualitatively, our findings are generally valid 
and representative of (unreconstructed) generic gold surfaces. The Maxwell relation, Eq. (O, 
is exact, and applies irrespective of the layer thickness. It rigorously shows that the sign of <j 
is negative, as long as the electron chemical potential grows with an increasing strain in the 
plane, e. This variation may be the superposition of several contributions a more detailed 
analysis of which has to be relegated to a forthcoming publication [30] . We mention here the 
most obvious terms, which cooperate to shift fi in the same direction, implying dfi/de>0. 

In order to understand the first mechanism, recall that the surface dipole describes how 
the electronic charge relaxes at a crystal surface after cleavage. This relaxation increases the 
work function W [16]. In the laterally compressed (e < 0) material, the increase in electron 
density and the normal expansion of the surface dipole (due to the transverse elastic response 
of the lattice) cooperate to increase the dipole strength, thereby increasing W and reducing 

Secondly, the trend dfj,/de>0 may also be understood in the framework of the microscopic 
picture proposed in Ref. [12], which considers the relaxation of the ion cores. Just as excess 
electronic charge {Sq < 0) on the surface of Au acts to displace the ion cores outwards relative 
to their Wigner-Seitz cell, resulting in an tendency for in-plane contraction (Sf > 0), the 
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converse effect is expected when the surface is compressed laterally (6e < 0): the ion cores 
will relax outward, closer to the surface of the electronic charge distribution. This trend may 
equivalently be viewed as the formation of an ionic dipole, on top of the electronic one. It 
again acts to increase the net dipole strength {S^ < 0). 

Thirdly, we note that the strain has also important consequences for the band structure, 
where it induces a band shift and a band flattening. For instance, the strain induced shift 
of the d-band center of gold relative to ep is large. It can exceed the shift of ^ found in our 
work by one or two orders of magnitude [6]. However, of interest to us is the change in ep 
with respect to the vacuum level. For this quantity, effects of band shift and flattening nearly 
compensate each other. 

The above effects depend on the atomic structure of the surface and, hence, on the ori- 
entation. This is consistent with the considerable difference in the values of ^ for (111) and 
(100). The observed trend, larger response on (111), agrees qualitatively with that of the 
surface stretch response to surface charge reported in Ref. [12], consistent with the suggested 
coupling between surface stress and stretch. 

The (110) surface presents a special case. In principle, it is possible that its parabolic 
n{e) carries information about real metallic films. However, we feel that our results arc not 
conclusive here because of the issues of convergence. It is conceivable that the slow convergence 
with film thickness signals the instability of the unreconstructed surface studied in our work 
with respect to a (missing row) reconstruction. An alternative reason could relate to the low 
symmetry of the (110) surface: the surface stress and its response to charging may well be 
anisotropic. By definition, the scalar is the mean of the two diagonal elements <;xx a-nd (^yy 
of a tensor response parameter. The small value of <r which might be extracted from the data 
in Fig. [ijd) could therefore be due to mutual cancelation of larger contributions <ixxSyy 

of 

opposite sign. 

Electrochemical experiments on Au (111) suggest values of <r in the range -0.67 to -0.9 
V [10,31]. The most negative values are observed for more weakly adsorbing electrolytes, 
where a smaller fraction of the total charge q is transferred away from the metal surface and 
into bonds with adsorbed ions [10]. In fact, experiments on Ft surfaces find ^ to decrease 
from -1.0 V to apparent saturation at -1.8 V when the adsorption is progressively reduced 
in dilution series [11]. The "intrinsic" magnitude of - in the absence of adsorption - may 
therefore be larger by as much as a factor of two as compared to experimental values recorded 
at typical concentrations. In this respect it is also remarkable that experimental data for the 
variation of the potential of zero charge of Au(lll) as a function of lateral strain - analogous 
to the present numerical procedure of computing d^/de - have been invoked to estimate the 
value of —qQ^dfi/de at -2 V [9], similar to our numerical result. We therefore conclude that 
our result for Au(lll) is consistent with the - limited - experimental data available at present. 
In view of the lack of fundamental understanding it is noteworthy that the computation for Au 
in vacuum reproduces the experimental sign of of Au in electrolytes, and that our results for 
<j are of the correct order of magnitude to explain the experimentally observed volume changes 
of nanoporous materials upon charging [1,2]. This supports the - not undisputed [7,9, 11] - 
notion of a large intrinsic value of related to processes and forces within the metal surface 
rather than to adsorption. 

In conclusion, an efficient way to compute the surface stress response to surface charg- 
ing by investigation of the electron chemical potential of strained but uncharged surfaces is 
presented and its application to Au surfaces by means of ab initio DFT with LDA is demon- 
strated, of Au (111) is evaluated to be -1.86 V, which is in reasonably good agreement with 
experimental observations. As compared to other capillary parameters, such as the surface 
energy, surface stress, and work function, the response is found to vary more strongly with 
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the surface orientation, the value of <; of the (100) surface being only about half of that of the 
(111) surface. 
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